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We analyze the spectral density of a single level quantum dot coupled to superconducting leads 
focusing on the Andreev states appearing within the superconducting gap. We use two complemen- 
tary approaches: the numerical renormalization group and the Hartree-Fock approximation. Our 
results show the existence of up to four bound states within the gap when the ground state is a spin 
doublet (-7T phase). Furthermore the results demonstrate the reliability of the mean field descrip- 
tion within this phase. This is understood from a complete correspondence that can be established 
between the exact and the mean field quasiparticle excitation spectrum within the gap. 
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INTRODUCTION 

Much progress has been achieved in recent years on 
the transport properties of quantum dots coupled to su- 
perconducting leads (for a review see [if). A central con- 
cept in the understanding of these properties is that of 
Andreev bound states (ABS), i.e. the bound states ap- 
pearing within the superconducting gap due to multiple 
Andreev reflections at the dot superconductor interfaces. 
When the two leads are superconducting the ABS depend 
on the superconducting phase difference and are thus cur- 
rent carrying states (typically most of the Josephson cur- 
rent is carried by the ABS PI). The interest in the ABS 
in these kind of systems has been increased by recent ex- 
periments allowing for their direct measurement through 
tunnel spectroscopy on carbon nanotube and graphene 
quantum dots @, Hj]. The issue is also related to the 
strong activity in the search of Majorana fermions, which 
would manifest as midgap states in different type of hy- 
brid nanostructures [J]. 

The simplest situation for analyzing the ABS spectrum 
is that of a single quantum dot (QD) with large energy 
level spacing which can be appropriately described by 
the single level Anderson model [B| . The presence of sub- 
gap states for the case of a magnetic impurity in a BCS 
superconducting host was already demonstrated in Refs. 
[S 0]- This model was afterward extended to analyze 
the Josephson transport properties and the so called 0- 
7T transition which signals the transition from a singlet 
(S = 0) to a doublet (S = 1/2) ground state HQ. In 
spite of these theoretical efforts the knowledge about the 
detailed structure of the ABS spectrum appears to be 
somewhat disperse in the literature. Thus, in a method 
like the non-crossing approximation (NCA), which is able 
to describe the 0-tt transition in the large charging en- 
ergy regime, a single subgap resonance appears whose 
crossing of the Fermi level signals the transition 15j . On 
the other hand, other approximations like Hartree-Fock 
(HFA) or exact diagonalizations in the infinite A limit, 



where A denotes the superconducting gap parameter in 
the leads, point to the existence of up to 4 levels symmet- 
rically located around the Fermi energy in the 7r-phase 
[lol \wL . The bound state spectrum has also been an- 
alyzed u sing the numerical renormalization group (NRG) 



method [18H2C§ although in these works only up to two 
ABS were identified. A previous NRG calculation show- 
ing up to four ABS exists [HJ although in this work only 
the single lead case was considered (i.e. without a phase 
difference between the leads). Taking into account all 
this rather fragmented evidence it seems worthwhile to 
investigate this issue in more detail using numerically ex- 
act results compared to different approximations. This 
is further motivated by the possibility of a direct experi- 
mental test along the lines of recent works [U, Q . 

In the present work we give a detailed analysis of the 
ABS for the single level Anderson model coupled to su- 
perconducting leads focusing in the regime Tr < A, 
where Tk denotes the Kondo temperature, which appears 
to be the relevant one for describing the experimental re- 
sults of Ref. Q ■ We use NRG calculations and compare 
the results with the mean-field approach provided by the 
HFA. It is found that when the system undergoes the 
transition to the 7r-phase there appear in general up to 4 
ABS in agreement with the analysis of the simple spin- 
polarized HFA [10] . Indeed, our analysis shows that the 
HFA provides a quite fair description of the bound state 
spectra except in the regime where Kondo correlations 
dominate over the superconducting ones. 

The rest of the paper is organized as follows: in Sect. 
□ we describe the model used for a QD coupled to su- 
perconducting leads and the basic theoretical analysis of 
its spectral properties; in Sect. „ we consider the simpler 
situation of a QD coupled to a single superconducting 
electrode (S-QD case) and present results for the ABS 
spectrum using both NRG and HFA, analyzing the range 
of validity of this last approximation for this case. In 
Sect. n this analysis is extended to a phase-biased S-QD- 
S system where we study in particular the behavior of 
the ABS spectrum around the transition between singlet 
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and doublet ground states. Finally in Sect, 
some concluding remarks. 



we give 



MODEL AND BASIC THEORETICAL ANALYSIS 



ferred to the lead chemical potential, i.e. = ek, u —^v) 
and A„ = | A„| exp (i<f> v ) is the (complex) superconduct- 
ing order parameter on lead v. Finally, Ht describes the 
coupling between the QD level to the leads and has the 
form 



A minimal model for a QD coupled to metallic elec- 
trodes in the regime where the energy level spacing Se is 
sufficiently large to restrict the analysis to a single spin- 
degenerate level is provided by the single level Anderson 
model with the Hamiltonian H = Ht + Hr + Ht + 
Hqd where Hqt> corresponds to the uncoupled dot given 
by 



Hqd — 



(1) 



where c' 0(7 creates and electron with spin a on the dot 
level located at eo and U is the local Coulomb interac- 
tion for two electrons with opposite spin within the dot 



(3) 



The coupling to the leads is usually characterized by a 
single parameter — 7rp„|V„| 2 , determining the width 
of the one-electron resonance. In this expression V v cor- 
responds to an average over the Fermi surface of Vk,v 
and p v denotes the corresponding density of states on 
the leads. 

In this work we are interested in the spectral properties 
which can be extracted from the dot Green's functions 
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On the other hand, H L . R describe the defined as G^, t') = -<*(*-*)<[*,(*). **(*)]+>. where 



(n 0a 

uncoupled left and right leads which are superconductors 
represented by a BCS Hamiltonian of the type 

H v = E&." C L,., < W + E { A » c il-.y c -H,v + h 

ka k 

(2) 

where c ka v creates an electron with spin a at the single- 
particle energy level of the lead v = L,R (usually re- 



= (c 0o 



r ) and the average is taken over the sys- 



tem ground state. In the case of the Anderson model with 
superconducting leads two types of ground states appear 
depending on the parameters: a non-degenerate S = 
ground state (0 phase) or a double-degenerate 5=1/2 
ground state (7r-phase). The Green's functions can be 
formally written in frequency space using the Lehmann 
representation 
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cj - (E m - E ) + i0+ 



{E m - E ) + 



(4) 
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where the system ground state, denoted by |$o,S*)) 
maybe degenerate (S z = ±1/2) and m labels the ex- 
cited states having an extra quasiparticle with respect 
to the ground state. In the degenerate case the total 
Green function is finally obtained as G r a = \ J2s g ■ 
The formal expression of Eq. ([J} allows a direct calcula- 
tion of the quasiparticle spectral densities using numer- 
ical methods like the NRG method which we describe 
further below. 

It would be interesting to compare the NRG results 
for the spectral density with those provided by different 
approximations specially with HFA which could provide 
a rather simple scheme for describing recent experimen- 
tal results on the ABS spectrum |2j. In this approxima- 



tion the dot Green's function is given by (Gp ) 1 = 
(G",(o))-i _ ±hf^ where Qr.,(o) is the non-interacting dot 

Green's function in Nambu space and the self-energy 
corresponds to the first order diagrams in the 
Coulomb interaction and are given by 



= ^22 F -a = U < «„,- 
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In the HFA both < no, a > and < co-fCoj, > have to be 
calculated self-consistently. The explicit expression for 

QT,HF is 
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hf = ( w - e - Tg(ui) - £^f ff 

^ rcos|/H-E^ 



r cos |/(c , 

e - rg(w) 
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^22,0- 
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where = Tfi = T/2, denotes the superconducting- 
phase difference and g(tu) = — A/(w)/u) = 
—uj/^/A 2 — u> 2 are the dimensioniess BCS Green's 
functions of the uncoupled leads. 

Within the HFA the transition from the to the 7r 
phase [1, [j| [Io[ is signaled by the existence of a spin 
broken symmetry solution with < n.o )(T >7^< "-o.-o- >• 
Although this symmetry breaking is not actually present 
in the exact solution, the HFA provides a very accurate 
description of the ABS spectrum within the n phase as 
will be shown along this work. This can be qualitatively 
understood by analyzing the low energy pole structure 
of the exact Green's function given by Eq. (jlj when the 
system is in the 7r-phase. In this case the two-fold degen- 
erate ground state can be labeled by S z = ±1/2. Quasi- 
particle excitations over this ground state can correspond 
to transitions to states with total spin either S = or 
5 = 1. However, in the last case the excitation energy is 
necessarily larger than A as these excited states involve 
an unpaired electron in the leads. Therefore the excita- 
tions within the gap can only arise from transitions to 
states with total spin equal to zero. It is then straight- 
forward to see that subgap electron-like excitations with 
spin a can only be created from the ground state with 
S z = —a while the hole-like excitations arise from the 
ground state S z = a. This structure is illustrated in 
Fig. n where we show schematically the subgap poles in 
G (7 ^s z ( UJ ) in the 7r-phase for S z = ±1/2. This is precisely 
the structure of the subgap excitations which are found 
in the spin-polarized HFA: the solutions for the major- 
ity and minority spin populations have only hole-like or 
electron-like character respectively. Therefore one can es- 
tablish a correspondence between the exact and the HFA 
excitations for the subgap states in the 7r phase. This 
correspondence based on the separation of electron and 
hole-like excitations is specially clear in the large A limit 
which we discuss in what follows. 

A — > oo case: As shown in previous works 13, 
the problem can be exactly diagonalized in this limit 
which already illustrates in a simple way the O-tt tran- 
sition. The states can be classified according to the 
total spin 5 = or 5 = 1/2. In the 5 = 1/2 sec- 
tor the energy levels are simply eo (doubly degener- 
ate) while in the 5 = case the states arc given by 
Eo,± = eo + U/2 ± y/(e Q + U/2) 2 + T 2 cos 2 <j>/2, leading 
to a phase transition for eo > Bo,-- Thus, in the 7r- 
phase the spectral density contains four ABS located at 
±U/2 ± ^(eo + U/2) 2 + F 2 cos 2 <j>/2. 

It is quite straightforward to see that this spectrum is 
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FIG. 1: (Color online) Subgap pole structure of the Green's 
functions components G a ,s z (w) in the 7r-phase for the two dif- 
ferent orientations of the ground state spin S z . As discussed 
in the text these excitations have either a purely hole-like or 
electron-like character, which allows a correspondence with 
the results of the HFA. Spin-symmetry is recovered when tak- 
ing the trace over the degenerate ground state. 



recovered exactly by the HFA. Indeed, in this limit the 
self-consistent HF solution is 
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r,HF 



eo + (2ff-l)¥ 



r cos § 



Tcosf 



eo 



(7) 



(2a + 1)% 

which exhibits the same spectrum as the exact solution. 
As in the 7r-phase U/2 > ^(eo + U/2) 2 + T 2 cos 2 (f>/2 the 
excitations for a = 1/2 have a hole- like character while 
those for a = —1/2 have an electron character. There- 
fore, it is not surprising that this approximation provides 
a rather good description of the ABS spectrum in the tt- 
phase for the full model. 

ABS SPECTRUM FOR THE S-QD CASE 

Before discussing the general case with two S leads 
and fixed phase difference it is worth analyzing the sim- 
pler case of an Anderson impurity coupled to a single 
BCS lead. This model exhibits also a transition to a 
degenerate 5=1/2 ground state and its spectral proper- 
ties are relevant to understand the transport properties 
in N-QD-S systems when < A [22]. To obtain the 



numerically exact ABS spectrum for this case we have 
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implemented an NRG algorithm following the lines of 
Refs. 



11, 21 



[231]. The idea behind the method is to 
discretize the energy levels in the leads on a logarithmic 
grid of energies A~" (with the dimensionless parameter 
A > 1 and 1 < n < N — > oo) with exponentially high 
resolution on the low-energy excitations. This discretiza- 
tion allows then to map the impurity model into a lin- 
ear "tight-binding" chain with hopping matrix elements 
decaying as A~™/ 2 with increasing site index n. The se- 
quence of Hamiltonians which is constructed by adding 
a new site in the chain is then diagonalized iteratively. 
As the number of states grows exponentially an adequate 
truncation scheme is required. 

The cutoff parameter A (as defined originally in Ref. 
[23 |) is chosen in order to ensure convergence of the spec- 
tral properties inside the gap. Depending on the value 
of U/T and the ratio Tk/A, we have chosen A varying 
between 2 and 4. For most of the results shown below we 
have checked that the value A = 4 provides already well 
converged results. In all the cases the usual correction 
Tnrg = A A T, where 



As 



1 1 + 1/A 
21- 1/A 



In A 



is used in order to correctly reproduce the exact A — > 1 
limit 2lJ, l23j, |24|. On the other hand, the maximum 
number of states N c kept in the iterative NRG procedure 
vary between ~ 300 in the S-QD case to ~ 800 for the 
S-QD-S case. 

In typical experiments the charging energy U adopts a 
nearly fixed value U > T while the dot level position can 
be varied. We thus first analyze the evolution of the ABS 
spectrum as a function of eo for fixed U /A and T/A. Fig. 
[2]shows the ABS spectrum and the corresponding weights 
for U/A = 1 and T/A = 0.2 obtained both within the 
NRG method and the HFA. The weights of the spectral 
density are calculated from the residues of the Green's 
functions (Eqs. dU) and ©) at the poles corresponding 
to the ABS energies. 

It is first worth noticing that the spectrum is character- 
ized by the presence of 4 ABS in the region \e + U/2\ < 
U/2, which corresponds to the S — 1/2 ground state, 
while only two ABS are present outside this region where 
the ground state is a singlet. As can be observed the 
HFA fairly reproduces not only the level positions but 
also their weights. Notice that the weights represented 
in the lower panels of Figs. [21 [3] correspond only to the 
electron like excitations which explains the asymmetry 
between positive and negative eo values. 

With increasing U the outermost ABS within the 
5=1/2 phase gradually approach the gap edge while its 
weight is reduced. Eventually, these states disappear for 
U/A ~ 2, as can be noticed in the inset of Fig. [5] which 
corresponds to the symmetric case. The ABS spectrum 




FIG. 2: (Color online) Upper panel: ABS spectrum of the 
S-QD system for the case U/A = 1 and T/A = 0.2 as a func- 
tion of the dot level position eo, calculated using the NRG 
(full lines) and HFA (dashed lines) respectively. The 7r-phase 
appearing around eo/A ~ —0.5 is characterized by the cross- 
ing of the internal Andreev levels and the presence of two 
extra states around E ~ ±0.5A, which disappear within the 
0-phase. The lower panel depicts the corresponding weights in 
the spectral density for the electron-like quasiparticle excita- 
tions. The inset shows the behavior of the ABS as a function 
of U/A in the electron-hole symmetric case eo = — U/2. The 
arrows in this plot indicate the U / A cases shown in the main 
frames of Figs. [2] and [3] respectively. 



properties for U/A = 4, illustrated in Fig. |31 clearly ex- 
hibits only two ABS within the gap. Again, as in the 
U /A = 1 case the agreement between the NRG results 
and the HFA is quite satisfactory. The main difference 
between both results appears at the crossing points be- 
tween the magnetic and non-magnetic regions where the 
HFA exhibits a small discontinuity. This discontinuity is 
due to the coexistence around these points of both types 
of solution in the HFA. In the results represented in Figs. 
[3] and [4] only the most stable HFA solution is shown. 

As a general remark one could state that the HFA 
reproduces fairly well the NRG results for arbitrary 
dot occupancy as far as the Kondo temperature Tk = 
^/[/r/2exp (— nU/8r) of the e-h symmetric case is 
smaller than A. Deviations with respect to the NRG re- 
sults could be expected when A/Tk becomes sufficiently 
small. This is illustrated in Fig. 2] where the ABS spec- 
trum is shown for the e-h symmetric case as a function of 
A/T K . As can be observed, when A/T K < 5 the HFA re- 
sults clearly deviates from the NRG ones as it predicts a 
magnetic solution up to the limit A/Tk — > (i.e. deep in 
the Kondo regime) whereas the NRG result becomes non- 
magnetic for A/Tk ~ 2.8. In contrast, in the opposite 
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FIG. 3: (Color online) Upper panel: ABS spectrum of the 
S-QD system for U/A = 4 and T/A = 0.2 as a function 
of the dot level position eo calculated using the NRG (full 
lines) and HFA (dashed lines) respectively. Notice that the 
outermost Andreev bound states have already merged with 
the continuum for this choice of parameters. The lower panel 
depicts the corresponding spectral weight for the electron-like 
excitations. 



limit A/Tk 3> 1 both results converge asymptotically to 
the A — > oo spectrum discussed in the previous section. 
We should point out that for the U/F ratio used in Fig. 
13] the "universal" limit (i.e. where all quantities depend 
only on the ratio A/Tk) is still not reached. For larger 
U/F ratios the transition occurs at smaller A/Tk, con- 
verging to a value ~ 1.7 when U/F > 10. This value is 
similar to the one reported in 2l| but somewhat larger 
than the one of Ref. [l2| obtained using quantum Monte 
Carlo techniques. One should notice also the difference 

II ED], which 



in the definition of Tk used in Refs. 
sponds to 0.4107 times the one used in the present work 
and also in Refs. 11, 12, IH 23 1 (the relation between 
the two definitions can be found in [25T]). 



FIG. 4: (Color online) ABS spectrum of the S-QD system 
as a function of A/Tk, where Tk is the Kondo temperature 
for the electron-hole symmetric case taking U/T = 5. Red 
and blue (full) lines correspond to the NRG calculation while 
black (dashed) lines are the HFA results. 



HFA is in this regime fairly good, as was already evident 
in Fig. [5] (which corresponds to the cf> — in this plot). 
When traversing the transition (middle panel of Fig. [5]) 
the agreement is less satisfactory due to the fact that 
the HFA result fully corresponds to the phase whereas 
within NRG the system is in a mixed 0' state (i.e. a 
mixed phase of character at 4> = and ir character at 
<f> = 7r with the absolute minimum energy corresponding 
to <j) = 0, see Ref. @). Finally, when the level position is 
sufficiently low both approaches predict a phase and the 
agreement in the ABS spectrum becomes progressively 
quite satisfactory again. This is the same trend which 
can be observed in Figs. f2]and|31 



CONCLUSIONS 



ABS SPECTRUM FOR THE S-QD-S CASE 

We analyze in this section the behavior of the ABS 
spectrum as a function of the phase difference for the 
S-QD-S case. The results shown in Fig. [5] correspond 
to the case U — A and F = 0.2 A, already analyzed in 
the previous section, for different values of the dot level 
illustrating the transition from the ir to the phase. The 
upper panel of Fig. [5] corresponds to the electron-hole 
symmetric case where the system exhibits 4 ABS inside 
the gap (notice that in the figure only the two electron- 
like states are shown). The agreement between NRG and 



In this work we have analyzed the subgap spectral den- 
sity of a single dot coupled to superconducting leads with 
the aim of clarifying some of the features of the ABS 
which appear to be controversial in the literature. By 
means of numerically exact NRG calculations we have 
shown that in general up to 4 ABS appear when the 
ground state becomes magnetic, i.e. in the 7r-phase. 
Within this phase the four states eventually reduce to 
only two for increasing U/A. Although the states are lo- 
cated symmetrically with respect to the Fermi level the 
electron-hole symmetry is in general broken. We have 
shown that this behavior is adequately reproduced by the 
HFA for a broad range of parameters, except very close to 
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FIG. 5: Evolution of the phase-dependent ABS spectrum 
when varying the dot level position for U / A = 1 and 
T/A = 0.2 calculated using the NRG method (full lines) 
and the HFA (dashed lines). From upper to lower panel 
e /A = -0.5,-0.7,-0.95,-1.2. Notice that only the elec- 
tron part of the spectrum is shown and that a different scale 
on the vertical axis is used for each panel in order to better 
visualize the details of the curves. 



the transition regions between the different phases. This 
approximation, however, is unable to reproduce the cor- 
rect behavior in the strong Kondo regime when Tk 3> A. 
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